function cels_invert_metropolis_pdf, x, mu=mu, sigma=sigma, pdf=pdf, positive=positive

  positive=intarr(n_elements(x))
  p=fltarr(n_elements(x))

  wh=where(pdf eq 'LOGNORMAL', count)
  if count ne 0 then begin
    p[wh]=(1./(x[wh]*sigma[wh]*sqrt(2.*!pi)))*exp( -1.*(alog(x[wh]) - mu[wh])^2 / 2./(sigma[wh]^2))
    positive[wh]=1
  endif

  wh=where(pdf eq 'EXPONENTIAL', count)
  if count ne 0 then begin
    p[wh]=mu[wh]*exp(-1.*mu[wh]*x[wh])
    positive[wh]=1
  endif

  wh=where(pdf eq 'GAUSSIAN', count)
  if count ne 0 then begin
    p[wh]=(1./(sigma[wh]*sqrt(2.*!pi)))*exp(-1.*(x[wh] - mu[wh])^2 / 2./(sigma[wh]^2))
    positive[wh]=0
  endif
  
  wh=where((x lt 0.) and positive, count)
  if count ne 0 then p[wh]=0.
  
  return, p
  
end

; docformat = 'rst'
;
;+
;
; :Purpose:
;   Solve for parameters(/state) vector using Metropolis-Hastings algorithm
;
; :Inputs:
;   sim_name: (string) Simulation name
;   
;   PDFs of each type of parameter (initial condition, non-local region, local region) can be specified
;   in {sim_name}.params.  The default PDFs are Gaussina for IC and NLR and Exponential for LR.
;   
;   Only Gaussian, Lognormal and Exponential PDFs are accepted at present.  Others can be added by changing
;   both routines in this file
;
; :Requires:
;   Combined sensitivity from *cels_sensitivity.pro*, stored in {input_directory}/{sim_name}/sensitivity.sav
;   
;   Measurements in {input_directory}/{sim_name}/measurements/measurements.sav
;   
;   State definitions in {input_directory}/{sim_name}/state.sav
;   
;   Region definitions in {input_directory}/{sim_name}/regions.sav
;
; :Keywords:
;   run_name: (string, input) will be added to the output file names to differentiate a particular Metropolis inversion from another using the same measurements/sensitivies, etc.
;
; :Outputs:
;   Histograms are stored in {input_directory}/{sim_name}/metropolis/metropolis{run_name}.sav
;   
;   Solutions (mean, percentiles, etc.) are stored in {input_directory}/{sim_name}/metropolis/metropolis{run_name}_solution.sav
;   
; :History:
; 	Written by: Matt Rigby, MIT, Aug 19, 2011
;
;-
Pro cels_invert_metropolis, sim_name, run_name=run_name

  sim_name='SF6_EA'

  compile_opt idl2
  do_metropolis=1

  if keyword_set(run_name) eq 0 then run_name=''
  StartY=cels_get_parameter(sim_name, 'STARTY')
  EndY=cels_get_parameter(sim_name, 'ENDY')
  Country_Track_params=cels_get_parameter(sim_name, 'Country_track')
  IC_uncertainty=cels_get_parameter(sim_name, 'IC_uncertainty')
  NLR_uncertainty=cels_get_parameter(sim_name, 'NLR_uncertainty')
  LR_uncertainty=cels_get_parameter(sim_name, 'LR_uncertainty')
  IC_pdf=cels_get_parameter(sim_name, 'IC_pdf')
  NLR_pdf=cels_get_parameter(sim_name, 'NLR_pdf')
  LR_pdf=cels_get_parameter(sim_name, 'LR_pdf')

  seed=systime(/seconds)


  restore, cels_filestr(/Input, sim_name + '/country_track.sav')
  restore, cels_filestr(/Input, sim_name + '/state.sav')
  restore, cels_filestr(/Input, sim_name + '/sensitivity.sav')

  if (n_elements(country_track_fraction[*, 0]) ne n_elements(x_ap)) or $
    (total(StrUpCase(Country_track_params[sort(country_track_params)]) ne country_track) gt 0.) then begin
    message, 'Run CELS_COUNTRIES_TO_TRACK'
  endif

;  nStations_all=!name.nStations_all
;  nStations=!name.nStations

  statesize=n_elements(x_ap)
  nMeasure=n_elements(y)
  nCountries=n_elements(country_track)

  ;Define STANDARD DEVIATIONS OF PRIOR PDFS
  ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

  x_ap_sd=(lr_uncertainty*x_ap)
  x_ap_sd[0:n_IC-1]=(IC_uncertainty*x_ap[0:n_IC-1])
  x_ap_sd[n_IC:N_NLRT-1]=(NLR_uncertainty*x_ap[n_IC:n_NLRT-1])

  ;Define which parts of x follow which PDF
  ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  
  x_pdf=strarr(StateSize)
  if IC_pdf ne !null then begin
    x_pdf[0:N_IC-1]=IC_PDF
  endif else begin
    print, 'CELS_INVERT_METROPOLIS, assuming GAUSSIAN IC'
    x_pdf[0:N_IC-1]='GAUSSIAN'
  endelse

  if NLR_pdf ne !null then begin
    x_pdf[N_IC:N_IC+N_NLRT-1]=NLR_PDF
  endif else begin
    print, 'CELS_INVERT_METROPOLIS, assuming GAUSSIAN NLR'
    x_pdf[N_IC:N_IC+N_NLRT-1]='GAUSSIAN'
  endelse

  if LR_pdf ne !null then begin
    x_pdf[N_IC+N_NLRT:-1]=LR_PDF
  endif else begin
    print, 'CELS_INVERT_METROPOLIS, assuming EXPONENTIAL LR'
    x_pdf[N_IC+N_NLRT:-1]='EXPONENTIAL'
  endelse

  if total(x_pdf eq '') then begin
    message, 'Need to define pdfs'
  endif


  ;Set PDF parameters
  wh=where(x_pdf eq 'LOGNORMAL', count)
  if count ne 0 then begin
    mu[wh]=alog(x_ap[wh]) - 0.5*alog(1. + x_ap_sd[wh]^2/x_ap[wh]^2) ;see wikipedia lognormal http://en.wikipedia.org/wiki/Log-normal_distribution
    sigma[wh]=sqrt(alog((x_ap_sd[wh]/x_ap[wh])^2 + 1.))
  endif
  wh=where(x_pdf eq 'EXPONENTIAL', count)
  if count ne 0 then begin
    mu=1./x_ap  ;exp distribution parameter mu*exp(-mu*x)
  endif
  wh=where((x_pdf ne 'EXPONENTIAL') and (x_pdf ne 'LOGNORMAL'), count)
  if count ne 0 then begin
    mu=x_ap 
    sigma=x_ap_sd
  endif
  

  if do_metropolis then begin
    
    ;Sample prior distribution
    ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
    
    burn_in=50000L
    nIt=100000L
    nIt_sub=1000L
    stepsize=1.5

    nBins=5000

    hist=fltarr(nBins, StateSize)
    locations=fltarr(nBins, StateSize)

    hist_country=fltarr(nBins, nCountries)
    locations_country=fltarr(nBins, nCountries)

    hist_global=fltarr(nBins)
    locations_global=fltarr(nBins)
    
    x=double(x_ap)

    n0=reform(H##x - y)
    n0T=total((n0/y_error)^2)

;    gridsize=x_ap_sd/10.

    p0=cels_invert_metropolis_pdf(x, mu=mu, sigma=sigma, pdf=x_pdf)

    x_it=fltarr(nIt_sub, StateSize)
    x_it[0, *]=x

    reject=0L
    it_sub=0L
    t0=systime(/seconds)

    print, 'Begin burn in'
    
    ;BEGIN ITERATIONS
    ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
    
    for it=0L, nIt+burn_in-1L do begin

      dx=randomn(seed, StateSize)*x_ap_sd*stepSize
;      if total(((x + dx) lt x_min) gt 0.) then begin
;        wh=where(x+dx lt x_min)
;        dx[wh]=x_min[wh] - x[wh]
;      endif

      dx_2d=dblarr(StateSize, nMeasure)
      for n=0, nMeasure-1 do dx_2d[*, n]=dx

      dy=H*dx_2d

      p1=cels_invert_metropolis_pdf(x+dx, mu=mu, sigma=sigma, pdf=x_pdf)

      for xi=0L, Statesize-1L do begin
        n1=reform(n0+dy[xi, *])
        n1T=total((n1/y_error)^2)
        pT=alog(p1[xi]/p0[xi])-0.5*(n1T - n0T)

        if alog(randomu(seed)) le pT then begin
          ;ACCEPT
          p0[xi]=p1[xi]
          x[xi]+=dx[xi]
          n0=n1
          n0T=n1T
        endif else begin
          ;REJECT
;          x[xi]=x0[xi]
          if it gt burn_in then reject+=1L
        endelse
      endfor
      

      if it ge burn_in then begin
        if ((it+1L) mod nIt_sub) eq 0 then begin

          for xi=0, StateSize-1 do begin
            hist[*, xi]+=histogram(x_it[*, xi], min=x_ap[xi] - 10.*x_ap_sd[xi], binsize=(50.*x_ap_sd[xi])/float(nBins), $
              nBins=nBins, locations=loctemp)
            locations[*, xi]=loctemp
          endfor

          ;NATIONAL EMISSIONS AVERAGED OVER ALL YEARS
          for xi=0, nCountries-1 do begin
            x_it_country=fltarr(nIt_sub)
            for n=0, nIt_sub-1 do x_it_country[n]=total(x_it[n, *]*country_track_fraction[*, xi])/(EndY - StartY)
            hist_country[*, xi]+=histogram(x_it_country, min=0., binsize=(total(x_ap[n_IC:-1]))/float(nBins), $
              nBins=nBins, locations=loctemp)
            locations_country[*, xi]=loctemp
          endfor

          x_it_global=fltarr(nIt_sub)
          for n=0, nIt_sub-1 do x_it_global[n]=total(x_it[n, n_IC:-1])/(EndY - StartY)
          hist_global+=histogram(x_it_global, min=0., binsize=(total(5.*x_ap[n_IC:-1]))/float(nBins), $
            nBins=nBins, locations=loctemp)
          locations_global=loctemp

          x_it[0, *]=x
          it_sub=0
          print, systime(/seconds) - t0, It - burn_in, 'acceptance: ', 1. - float(reject)/float(It-burn_in+1L)/float(StateSize)
        endif else begin
          x_it[it_sub+1, *]=x
          it_sub+=1L
        endelse
      endif else begin
        x_it[0, *]=x
        reject=0L
      endelse
      
      
    endfor

    save, filename=cels_filestr(/Input, sim_name + '/metropolis/metropolis' + run_name + '.sav'), hist, locations, $
      nIt, reject, burn_in, x_ap, x_ap_sd, hist_global, locations_global, hist_country, locations_country

  endif else restore, cels_filestr(/Input, sim_name + '/metropolis/metropolis' + run_name + '.sav')

  print, 'Acceptance ratio: ', 1. - float(reject)/float(nIt)/float(StateSize)

  x_mean=fltarr(StateSize)
  x_ml=fltarr(StateSize)
  x_sd=fltarr(StateSize)
  x_50=fltarr(StateSize)
  x_16=fltarr(StateSize)
  x_84=fltarr(StateSize)
  x_25=fltarr(StateSize)
  x_75=fltarr(StateSize)

  Country_mean=fltarr(nCountries)
  Country_ml=fltarr(nCountries)
  Country_sd=fltarr(nCountries)
  Country_50=fltarr(nCountries)
  Country_16=fltarr(nCountries)
  Country_84=fltarr(nCountries)
  Country_25=fltarr(nCountries)
  Country_75=fltarr(nCountries)

  
  ;Calculate mean, percentiles, max_likelihood
  nBins=n_elements(hist[*, 0])
  
  nIt=total(hist[*, 0])
  
  ;Regional
  ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  for xi=0, StateSize-1 do begin
    dlocations=locations[1, xi] - locations[0, xi]
    x_mean[xi]=total(hist[*, xi]*(locations[*, xi] + dlocations/2.))/total(hist[*, xi])
    x_sd[xi]=sqrt(total(hist[*, xi]*(locations[*, xi] + dlocations/2. - x_mean[xi])^2)/total(hist[*, xi]))
    cum=0.
    nIt_hist=total(hist[*, xi])
    for n=0L, nBins-1 do begin
      cum+=hist[n, xi]
      if cum lt nIt_hist then begin 
        if cum gt 0.16*float(nIt) and x_16[xi] eq 0. then x_16[xi]=locations[n, xi]; + dlocations/2.
        if cum gt 0.5*float(nIt) and x_50[xi] eq 0. then x_50[xi]=locations[n+1, xi]; + dlocations/2.
        if cum gt 0.84*float(nIt) and x_84[xi] eq 0. then x_84[xi]=locations[n+1, xi]; + dlocations/2.
        if cum gt 0.25*float(nIt) and x_25[xi] eq 0. then x_25[xi]=locations[n+1, xi]; + dlocations/2.
        if cum gt 0.75*float(nIt) and x_75[xi] eq 0. then x_75[xi]=locations[n+1, xi]; + dlocations/2.
      endif
    endfor
    dummy=max(hist[*, xi], maxI)
    x_ml[xi]=locations[maxI, xi] + dlocations/2.
  endfor

  ;National
  ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

  for xi=0, nCountries-1 do begin
    dlocations_country=locations_country[1, xi] - locations_country[0, xi]
    Country_mean[xi]=total(hist_country[*, xi]*(locations_country[*, xi] + dlocations_country/2.))/total(hist_country[*, xi])
    Country_sd[xi]=sqrt(total(hist_country[*, xi]*(locations_country[*, xi] + dlocations_country/2. - Country_mean[xi])^2)/total(hist_country[*, xi]))
    cum=0.
    nIt_hist_country=total(hist_country[*, xi])
    for n=0L, nBins-1 do begin
      cum+=hist_country[n, xi]
      if cum lt nIt_hist_country then begin 
        if cum gt 0.16*float(nIt) and Country_16[xi] eq 0. then Country_16[xi]=locations_country[n, xi]; + dlocations_country/2.
        if cum gt 0.5*float(nIt) and Country_50[xi] eq 0. then Country_50[xi]=locations_country[n+1, xi]; + dlocations_country/2.
        if cum gt 0.84*float(nIt) and Country_84[xi] eq 0. then Country_84[xi]=locations_country[n+1, xi]; + dlocations_country/2.
        if cum gt 0.25*float(nIt) and Country_25[xi] eq 0. then Country_25[xi]=locations_country[n+1, xi]; + dlocations_country/2.
        if cum gt 0.75*float(nIt) and Country_75[xi] eq 0. then Country_75[xi]=locations_country[n+1, xi]; + dlocations_country/2.
      endif
    endfor
    dummy=max(hist_country[*, xi], maxI)
    Country_ml[xi]=locations_country[maxI, xi] + dlocations_country/2.
  endfor

  ;GLOBAL
  ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

  dlocations_global=locations_global[1] - locations_global[0]
  global_mean=total(hist_global*(locations_global + dlocations_global/2.))/total(hist_global)
  global_sd=sqrt(total(hist_global*(locations_global + dlocations_global/2. - global_mean)^2)/total(hist_global))
  cum=0.
  nIt_hist_global=total(hist_global)
  for n=0L, nBins-1 do begin
    cum+=hist_global[n]
    if cum lt nIt_hist_global then begin 
      if cum gt 0.16*float(nIt) and n_elements(global_16) eq 0 then global_16=locations_global[n]; + dlocations_global/2.
      if cum gt 0.5*float(nIt) and n_elements(global_50) eq 0 then global_50=locations_global[n+1]; + dlocations_global/2.
      if cum gt 0.84*float(nIt) and n_elements(global_84) eq 0 then global_84=locations_global[n+1]; + dlocations_global/2.
      if cum gt 0.25*float(nIt) and n_elements(global_25) eq 0 then global_25=locations_global[n+1]; + dlocations_global/2.
      if cum gt 0.75*float(nIt) and n_elements(global_75) eq 0 then global_75=locations_global[n+1]; + dlocations_global/2.
    endif
  endfor
  dummy=max(hist_global, maxI)
  global_ml=locations_global[maxI] + dlocations_global/2.

  save, filename=cels_filestr(/Input, sim_name + '/metropolis/metropolis' + run_name + '_solution.sav'), $
    hist, locations, nIt, reject, x_mean, x_ml, x_50, x_16, x_84, x_sd, x_ap, x_ap_sd, x_25, x_75, $
    Global_mean, Global_ml, Global_50, Global_16, Global_84, Global_sd, Global_25, Global_75, $
    Country_mean, Country_ml, Country_50, Country_16, Country_84, Country_sd, Country_25, Country_75

end
